In [ ]:
import warnings
import pandas as pd
import numpy as np
import random
import seaborn as sns
import matplotlib.pyplot as plt

from os.path import join
from IPython import display
from sklearn.datasets import load_digits
from sklearn.datasets import make_blobs
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.metrics import silhouette_score, homogeneity_score, completeness_score, v_measure_score # и другие метрики
from sklearn.metrics.cluster import adjusted_rand_score, adjusted_mutual_info_score
from sklearn.cluster import KMeans # а также другие алгоритмы
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics.pairwise import chi2_kernel
from sklearn.cluster import Birch, DBSCAN, BisectingKMeans, MiniBatchKMeans, SpectralClustering
In [ ]:
DATA_PATH = "../data"
plt.rcParams["figure.figsize"] = 12, 9
sns.set_style("whitegrid")
warnings.filterwarnings("ignore")

SEED = 111
random.seed(SEED)
np.random.seed(SEED)

Задание 1. Реализация Kmeans¶

5 баллов

В данном задании вам предстоит дописать код класса MyKMeans. Мы на простом примере увидим, как подбираются центры кластеров и научимся их визуализировать.

Сгенерируем простой набор данных, 400 объектов и 2 признака (чтобы все быстро работало и можно было легко нарисовать):

In [ ]:
X, true_labels = make_blobs(400, 2, centers=[[0, 0], [-4, 0], [3.5, 3.5], [3.5, -2.0]])

Напишем функцию visualize_clusters, которая по данным и меткам кластеров будет рисовать их и разукрашивать:

In [ ]:
def visualize_clusters(X, labels):
    """
    Функция для визуализации кластеров
        :param X: таблица объекты х признаки
        :param labels: np.array[n_samples] - номера кластеров
    """
    
    unique_labels = np.sort(np.unique(labels))
    sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=labels, 
                    palette=sns.color_palette("RdPu")[1:], legend=False,
                    hue_order=unique_labels)
    plt.xlabel("$X_1$", fontsize=18)
    plt.ylabel("$X_2$", fontsize=18)
    
    for label in labels:
        center = X[(labels == label)].mean(axis=0)
        plt.scatter(center[0], center[1], s=80, c="#201F12", marker=(5, 2))
In [ ]:
visualize_clusters(X, true_labels)

Напишем свой класс MyKMeans, который будет реализовывать алгоритм кластеризации K-средних. Напомним сам алгоритм:

  1. Выбераем число кластеров (K)
  2. Случайно инициализируем K точек (или выбираем из данных), это будут начальные центры наших кластеров
  3. Далее для каждого объекта считаем расстояние до всех кластеров и присваиваем ему метку ближайщего
  4. Далее для каждого кластера считаем "центр масс" (среднее значение для каждого признака по всем объектам кластера)
  5. Этот "центр масс" становится новым центром кластера
  6. Повторяем п.3, 4, 5 заданное число итераций или до сходимости

Во время предсказания алгоритм просто находит ближайщий центроид (центр кластера) для тестового объекта и возвращает его номер.

Реализуйте методы:

  • _calculate_distance(X, centroid) - вычисляет Евклидово расстояние от всех объектов в Х до заданного центра кластера (centroid)
  • predict(X) - для каждого элемента из X возвращает номер кластера, к которому относится данный элемент
In [ ]:
class MyKMeans:
    def __init__(self, n_clusters, init="random", max_iter=300, visualize=False):
        """
        Конструктор класса MyKMeans
            :param n_clusters: число кластеров
            :param init: способ инициализации центров кластеров
                'random' - генерирует координаты случайно из нормального распределения
                'sample' - выбирает центроиды случайно из объектов выборки
            :param max_iter: заданное число итераций 
                (мы не будем реализовывать другой критерий остановки)
            :param visualize: рисовать ли кластеры и их центроиды в процессе работы
                код будет работать сильно дольше, но красиво...
        """
        
        assert init in ["random", "sample"], f"Неизвестный метод инициализации {init}"
        self.n_clusters = n_clusters
        self.init = init
        self.max_iter = max_iter
        self.centroids = None
        self.visualize = visualize
       
    
    def fit(self, X):
        """
        Подбирает оптимальные центры кластеров
            :param X: наши данные (n_samples, n_features)
        :return self: все как в sklearn
        """
        
        n_samples, n_features = X.shape
        
        # Инициализация центров кластеров
        if self.init == "random":
            centroids = np.random.randn(self.n_clusters, n_features)
        elif self.init == "sample":
            centroids_idx = np.random.choice(np.arange(n_samples), 
                                             size=self.n_clusters, 
                                             replace=False)
            centroids = X[centroids_idx]
        
        # Итеративно двигаем центры
        for _ in range(self.max_iter):
            # Посчитаем расстояния для всех объектов до каждого центроида
            dists = []
            for centroid in centroids:
                dists.append(self._calculate_distance(X, centroid))
            dists = np.concatenate(dists, axis=1)
            # Для каждого объекта найдем, к какому центроиду он ближе
            cluster_labels = np.argmin(dists, axis=1)
            
            # Пересчитаем центр масс для каждого кластера
            centroids = []
            for label in np.sort(np.unique(cluster_labels)):
                center = X[(cluster_labels == label)].mean(axis=0)
                centroids.append(center)
            
            # Отрисуем точки, покрасим по меткам кластера, а также изобразим центроиды
            if self.visualize:
                visualize_clusters(X, cluster_labels)
                display.clear_output(wait=True)
                display.display(plt.gcf())
                plt.close()
                
        self.centroids = np.array(centroids)
        
        return self
    
    
    def predict(self, X):
        """
        Для каждого X возвращает номер кластера, к которому он относится
            :param X: наши данные (n_samples, n_features)
        :return cluster_labels: метки кластеров
        """
        
        """
        YOUR CODE IS HERE
        """
        dists = []
        for centroid in self.centroids:
            dists.append(self._calculate_distance(X, centroid))
        dists = np.concatenate(dists, axis=1)
        cluster_labels = np.argmin(dists, axis=1)
            
        return cluster_labels
        
        
    def _calculate_distance(self, X, centroid):
        """
        Вычисляет Евклидово расстояние от всех объектов в Х до заданного центра кластера (centroid)
            :param X: наши данные (n_samples, n_features)
            :param centroid: координаты центра кластера
        :return dist: расстояния от всех X до центра кластера
        """
        
        """
        YOUR CODE IS HERE
        """
        dist = np.sqrt(np.sum(np.square(X - centroid), axis=1)).reshape(-1, 1)
        
        return dist
    
    
    def __repr__(self):
        return f"Привет, я твой KMeans (/¯◡ ‿ ◡)/¯☆*"

Обучите MyKMeans на наших игручешных данных, добейтесь сходимости. Не забудьте поставить visualize=True, чтобы посмотреть на красивые картинки. Также попробуйте различные способы инициализации центроидов и скажите, какой лучше подошел в этой ситуации.

In [ ]:
initialization = "random"
iterations = 10

my_k_means_random = MyKMeans(n_clusters=4, visualize=True, init=initialization, max_iter=iterations)
my_k_means_random.fit(X)
pred_random = my_k_means_random.predict(X)

print(
    f'{initialization} initialization, {iterations} iterations:'
    '\n'
    f'The mean Silhouette Coefficient: {silhouette_score(X, pred_random)}'
    '\n'
    f'Rand index adjusted for chance: {adjusted_rand_score(true_labels, pred_random)}'
    '\n'
    f'Adjusted Mutual Information: {adjusted_mutual_info_score(true_labels, pred_random)}'
    )
random initialization, 10 iterations:
The mean Silhouette Coefficient: 0.5643117184043004
Rand index adjusted for chance: 0.9278697796079957
Adjusted Mutual Information: 0.9051395864810662
In [ ]:
initialization = "sample"
iterations = 10

my_k_means_sample = MyKMeans(n_clusters=4, visualize=True, init=initialization, max_iter=iterations)
my_k_means_sample.fit(X)
pred_sample = my_k_means_sample.predict(X)

print(
    f'{initialization} initialization, {iterations} iterations:'
    '\n'
    f'The mean Silhouette Coefficient: {silhouette_score(X, pred_sample)}'
    '\n'
    f'Rand index adjusted for chance: {adjusted_rand_score(true_labels, pred_sample)}'
    '\n'
    f'Adjusted Mutual Information: {adjusted_mutual_info_score(true_labels, pred_sample)}'
    )
sample initialization, 10 iterations:
The mean Silhouette Coefficient: 0.5643117184043004
Rand index adjusted for chance: 0.9278697796079957
Adjusted Mutual Information: 0.9051395864810661

Here I took intentionally low iteration number (10) to show that on these data all three metrics are slightly higher if clustaers were initialized from sample points. However, the difference is extremely small. Let's try on 50 iterations.

In [ ]:
initialization = "random"
iterations = 50

my_k_means_random = MyKMeans(n_clusters=4, visualize=True, init=initialization, max_iter=iterations)
my_k_means_random.fit(X)
pred_random = my_k_means_random.predict(X)

print(
    f'{initialization} initialization, {iterations} iterations:'
    '\n'
    f'The mean Silhouette Coefficient: {silhouette_score(X, pred_random)}'
    '\n'
    f'Rand index adjusted for chance: {adjusted_rand_score(true_labels, pred_random)}'
    '\n'
    f'Adjusted Mutual Information: {adjusted_mutual_info_score(true_labels, pred_random)}'
    )
random initialization, 50 iterations:
The mean Silhouette Coefficient: 0.5644544929974866
Rand index adjusted for chance: 0.9340183005202644
Adjusted Mutual Information: 0.9117474366398979
In [ ]:
initialization = "sample"
iterations = 50

my_k_means_sample = MyKMeans(n_clusters=4, visualize=True, init=initialization, max_iter=iterations)
my_k_means_sample.fit(X)
pred_sample = my_k_means_sample.predict(X)

print(
    f'{initialization} initialization, {iterations} iterations:'
    '\n'
    f'The mean Silhouette Coefficient: {silhouette_score(X, pred_sample)}'
    '\n'
    f'Rand index adjusted for chance: {adjusted_rand_score(true_labels, pred_sample)}'
    '\n'
    f'Adjusted Mutual Information: {adjusted_mutual_info_score(true_labels, pred_sample)}'
    )
sample initialization, 50 iterations:
The mean Silhouette Coefficient: 0.5644544929974866
Rand index adjusted for chance: 0.9340183005202644
Adjusted Mutual Information: 0.9117474366398978

As we can see, the more iterations, the less initialization type does affect the model quality.

Задание 2. Подбираем лучшую иерархическую кластеризацию¶

5 баллов

На лекции были рассмотрены различные расстояния, которые могут служить метриками различия между объектами. Также мы разобрали несколько алгоритмов кластеризации, в том числе и иерархическую. Часто может быть непонятно, какой алгоритм и какую метрику расстояния нужно взять. Давайте упростим себе жизнь и напишем функцию algorithm_selection, которая будет на основании переданных ей:

  • метрик расстояния (можно брать все, что было на лекциях, минимум 4)
  • параметра linkage ('average', 'single', 'complete')
  • и метрик качества кластеризации ('Homogeneity', 'Completeness', 'V-measure', 'Silhouette')

будет выводить итоговую таблицу, в которую войдут столбцы:

  • distance (метрика расстояния)
  • linkage (алгоритм иерархической кластеризации)
  • homogenity
  • completeness
  • v_measure
  • silhouette

В результате по этой таблице, нужно сделать вывод о том, какой алгоритм кластеризации и с какими гиперпараметрами лучше всего сработал.

Загрузим наши данные:

In [ ]:
data = load_digits()
X, y = data.data, data.target

Let's take a look at the shapes of X and y, how many classes are present and visualize them and see the representability of classes.

In [ ]:
y.dtype
Out[ ]:
dtype('int64')
In [ ]:
print(f'X shape: {X.shape}, y shape: {y.shape}, number of classes: {np.unique(y).shape[0]}, class labels: {np.unique(y)}')
X shape: (1797, 64), y shape: (1797,), number of classes: 10, class labels: [0 1 2 3 4 5 6 7 8 9]
In [ ]:
plt.imshow(X[0].reshape(8, 8).astype(int), cmap="RdPu")
plt.axis("off");
In [ ]:
fig, axes = plt.subplots(2, 5, figsize=(13.5, 6))

fig.suptitle("Representation of each class object", size=16, weight="bold", y=0.91)

class_nums = list(map(int, range(10)))

object_class = 0
for row in axes:
    for col in row:
        object_index = np.where(y == object_class)[0][0]
        image = X[object_index].reshape(8, 8)
        col.set_title(f'Class_{object_class}', weight='bold', size=12)
        col.set_axis_off()
        col.imshow(image, cmap='RdPu')
        object_class += 1
plt.subplots_adjust(wspace=-0.3)
fig.tight_layout(pad=2.0)
In [ ]:
_, object_count = np.unique(y, return_counts=True)

plt.figure(figsize=(15, 6))
sns.barplot(x=list(map(lambda x: f'Class_{class_nums.index(x)}', class_nums)), 
            y=object_count,
            color=sns.color_palette("RdPu")[3],
            edgecolor=sns.color_palette("RdPu")[4],
            alpha=0.85)
plt.xlabel("Class number and name", size=13, weight="bold")
plt.ylabel("Object count", size=13, weight="bold")
plt.title("Class representation", size=15, weight="bold")
plt.axhline(y=object_count.min(), color=sns.color_palette("RdPu")[5], linestyle='-.', linewidth=0.9)
plt.axhline(y=object_count.max(), color=sns.color_palette("RdPu")[5], linestyle='-.', linewidth=0.9)
plt.axhline(y=object_count.mean(), color=sns.color_palette("RdPu")[5], linestyle='-.', linewidth=1.2)
;
Out[ ]:
''

Работать будем с изображениями рукописных цифр. Тут все похоже на данные для 1 домашнего задания, каждая картинка представлена вектором из 64 элементов (числа от 0 до 255). Чтобы ее нарисовать мы должны сделать reshape в картинку 8 на 8 пикселей. Вам нужно будет выбрать наилучший способ кластеризации при помощи функции algorithm_selection, которую вы реализуете. Для некоторых метрик кластеризации требуются метки классов объектов (они хранятся в переменной y).

YOUR TURN TO CODE¶
In [ ]:
def algorithm_selection(X, distances, algorithms, metrics, cluster_number):
    """
    Для заданных алгоримов кластеризации и гиперпараметров 
    считает различные метрики кластеризации
        :param X: наши данные (n_samples, n_features)
        :param distances: список возможных метрик расстояния
        :param algorithm: параметр linkage ('average', 'single', 'complete')
        :param metrics: список возможных метрик качества кластеризации
        :param cluster_number: количество кластеров - добавила, чтобы потом для клеток использовать, и вообще норм было варьировать вводные данные
    :return compare_dataframe: таблица с метриками кластеризации
    """

    """
    YOUR CODE IS HERE
    """

    col_names = list(("distance", "linkage"))
    col_names.extend(sorted(metrics))
    compare_dataframe = pd.DataFrame()

    for distance in distances:
        for algorithm in algorithms:
            metric_result = list()

            if distance == "chi-square":
                chi2_distances = chi2_kernel(X)
                np.fill_diagonal(chi2_distances, 0)
                alg_clust = AgglomerativeClustering(
                    n_clusters=cluster_number, 
                    affinity='precomputed', 
                    linkage=algorithm)
                prediction = alg_clust.fit_predict(chi2_distances)
                if "silhouette" in metrics:
                    silhouette = silhouette_score(chi2_distances, prediction, metric='precomputed')
            else:

                alg_clust = AgglomerativeClustering(
                    n_clusters=cluster_number, 
                    affinity=distance, linkage=algorithm)
                prediction = alg_clust.fit_predict(X)
                if "silhouette" in metrics:
                    silhouette = silhouette_score(X, prediction, metric=distance)
            
            if "completeness" in metrics:
                completeness = completeness_score(y, prediction)
                metric_result.append(completeness)
            if "homogenity" in metrics:
                homogeneity = homogeneity_score(y, prediction)
                metric_result.append(homogeneity)
            if "silhouette" in metrics:
                metric_result.append(silhouette)
            if "v_measure" in metrics:
                v_measure = v_measure_score(y, prediction)
                metric_result.append(v_measure)

            new_row = list((distance, algorithm))
            new_row.extend(metric_result)
            compare_dataframe = compare_dataframe.append(pd.Series(new_row), ignore_index=True)
    compare_dataframe.set_axis(col_names, axis=1, inplace=True)
    
    return compare_dataframe
In [ ]:
distances = ["euclidean", "manhattan", "canberra", "chi-square", "cosine"]
linkages = ["average", "single", "complete"]
metrics = ["homogenity", "completeness", "v_measure", "silhouette"]
In [ ]:
algos_selection = algorithm_selection(X, distances, linkages, metrics, np.unique(y).shape[0])
algos_selection
Out[ ]:
distance linkage completeness homogenity silhouette v_measure
0 euclidean average 0.782342 0.655393 0.146625 0.713263
1 euclidean single 0.273755 0.005056 -0.133621 0.009929
2 euclidean complete 0.637772 0.590549 0.119271 0.613253
3 manhattan average 0.778942 0.566751 0.105238 0.656117
4 manhattan single 0.276237 0.005102 -0.203816 0.010019
5 manhattan complete 0.686109 0.646372 0.160617 0.665648
6 canberra average 0.779789 0.573860 0.108401 0.661161
7 canberra single 0.280506 0.006144 -0.157872 0.012025
8 canberra complete 0.561588 0.528702 0.104411 0.544649
9 chi-square average 0.003473 0.003367 -0.750211 0.003419
10 chi-square single 0.276560 0.005108 -0.994992 0.010031
11 chi-square complete 0.003609 0.003565 -0.721485 0.003587
12 cosine average 0.808762 0.648268 0.218229 0.719676
13 cosine single 0.273389 0.005050 -0.281789 0.009916
14 cosine complete 0.577415 0.528727 0.144179 0.551999
In [ ]:
algos_selection[["homogenity", "completeness", "v_measure", "silhouette"]].idxmax()
Out[ ]:
homogenity       0
completeness    12
v_measure       12
silhouette      12
dtype: int64
In [ ]:
algos_selection.loc[[12]]
Out[ ]:
distance linkage completeness homogenity silhouette v_measure
12 cosine average 0.808762 0.648268 0.218229 0.719676
In [ ]:
algos_selection.loc[[0]]
Out[ ]:
distance linkage completeness homogenity silhouette v_measure
0 euclidean average 0.782342 0.655393 0.146625 0.713263

Okay, so we see that most of the metrics except homogeneity score are best with following parameters: cosine distances + average linkage (homogeneity score is still not bad), except for homogeneity score which is better if euclidean distances + average linkage (completeness, silhouette and v_measure are still not bad).

So, I'd choose the hyperparameter combination between these two, preferrably, cosine distances + average linkage.

Задание 3. Аннотация клеточных типов¶

суммарно 10 баллов

В этом задании вам предстоит применять полученные знания о кластеризации для аннотации данных, полученных при помощи проточной цитометрии. Каждая строка в данных это одна клетка, столбцы FSC-A, SSC-A косвенным образом свидетельствуют о размере клеток, остальные показывают интенсивность флуоресценции различных поверхностных маркеров. Ваша задача определить для каждой клетки, к какому типу она относится.

3.1. EDA¶

1.5 балла

Проведите EDA:

  1. Посмотрите на данные (можно попробовать метод describe)
  2. Сколько у нас клеток / признаков
  3. Постройте распределения каждого признака (есть ли очевинные выбросы?)
In [ ]:
fc_data = pd.read_csv(join(DATA_PATH, "flow_c_data.csv"), index_col=0)
YOUR TURN TO CODE:¶
In [ ]:
fc_data.head()
Out[ ]:
FSC-A SSC-A HLA-DR CD3 CD45 CD14 CD56 CD19 CD11c
0 0.722199 0.154232 0.494449 0.512853 4.057328 1.019255 4.007273 -1.406912 0.384406
1 0.751073 0.114618 4.671924 0.822042 4.933149 0.159388 -0.216407 6.034845 0.577039
2 1.108188 0.150011 -0.422857 4.006353 5.616384 0.519223 0.783077 -0.318668 -0.493832
3 0.553310 0.119929 0.345968 3.332310 5.260935 1.481229 0.264846 0.694883 3.416960
4 0.774884 0.076984 4.984963 0.885501 4.562367 0.217117 0.086279 5.563345 0.471591
In [ ]:
fc_data.shape # 10000 cells, 9 features
Out[ ]:
(10000, 9)
In [ ]:
fc_data.isna().sum() # great, no NAs to care about
Out[ ]:
FSC-A     0
SSC-A     0
HLA-DR    0
CD3       0
CD45      0
CD14      0
CD56      0
CD19      0
CD11c     0
dtype: int64
In [ ]:
fc_data.describe()
Out[ ]:
FSC-A SSC-A HLA-DR CD3 CD45 CD14 CD56 CD19 CD11c
count 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000 10000.000000
mean 0.962458 0.165415 1.791579 3.280391 5.169845 1.097182 0.412223 1.219683 0.569898
std 0.500012 0.190773 2.027281 2.112079 0.518247 1.375876 1.176967 2.285305 1.502445
min 0.144786 0.021657 -7.634172 -2.517639 -8.329666 -8.833883 -12.358701 -9.640775 -11.085437
25% 0.798804 0.105550 0.275056 1.401943 4.922227 0.389746 -0.195978 -0.146948 -0.498781
50% 0.908574 0.130543 1.194313 4.373660 5.234119 0.672853 0.053267 0.244121 0.272297
75% 1.052034 0.162536 3.696976 4.842212 5.481218 1.065076 0.378576 0.908374 1.220576
max 42.954485 15.460997 12.946668 13.844391 13.710922 6.294236 6.461899 7.208721 6.840864

First, let's see how the data is distributed trying to identify outliers.

In [ ]:
fc_data.columns
Out[ ]:
Index(['FSC-A', 'SSC-A', 'HLA-DR', 'CD3', 'CD45', 'CD14', 'CD56', 'CD19',
       'CD11c'],
      dtype='object')
In [ ]:
fig, axes = plt.subplots(6, 3, figsize=(15, 15), gridspec_kw={"height_ratios": [1, 4] * 3})

idx=0
for i in list(range(6))[::2]:
    for j in range(3):
        axes[i][j].set_title(f'{fc_data.columns[idx]}', size=15, weight='bold')
        sns.boxplot(fc_data, x=fc_data.columns[idx], width=0.5, ax=axes[i][j], color=sns.color_palette("RdPu")[4], fliersize=1.5)
        sns.distplot(fc_data[fc_data.columns[idx]], ax=axes[i + 1][j], color=sns.color_palette("RdPu")[4])
        idx += 1
fig.subplots_adjust(hspace=0.9)

Some obviously outlying observations exclusion.

In [ ]:
outliers = fc_data[
    (fc_data["FSC-A"] > 30) 
    | (fc_data["SSC-A"] > 5)
    | (fc_data["HLA-DR"] < -5)
    | (fc_data["HLA-DR"] > 10)
    | (fc_data["CD3"] > 10)
    | (fc_data["CD45"] < 0)
    | (fc_data["CD45"] > 10)
    | (fc_data["CD14"] < -1)
    | (fc_data["CD56"] < -5)
    | (fc_data["CD19"] < -5)
    | (fc_data["CD11c"] < -3)
    ]
print(f'{outliers.shape[0]} outliers identified')
fc_data_cleared = fc_data.drop(outliers.index)
11 outliers identified

Have to see feature distribution again. Looks better in terms of outliers, obviously not normally distributed though but anyway it's not crucial for clustering.

In [ ]:
fig, axes = plt.subplots(6, 3, figsize=(15, 15), gridspec_kw={"height_ratios": [1, 4] * 3})

idx=0
for i in list(range(6))[::2]:
    for j in range(3):
        axes[i][j].set_title(f'{fc_data_cleared.columns[idx]}', size=15, weight='bold')
        sns.boxplot(fc_data_cleared, x=fc_data_cleared.columns[idx], width=0.5, ax=axes[i][j], color=sns.color_palette("RdPu")[4], fliersize=1.5)
        sns.distplot(fc_data_cleared[fc_data_cleared.columns[idx]], ax=axes[i + 1][j], color=sns.color_palette("RdPu")[4])
        idx += 1
fig.subplots_adjust(hspace=0.9)

3.2. Кластеризация¶

4.5 балла

При ручной аннотации клеточных типов обычно поступают следующим образом:

  1. При помощи методов понижения размерности рисуют наши наблюдения, чтобы примерно оценить число клеточных типов
  2. Проводят кластеризацию наблюдений (для некоторых методов нужно заранее задать число кластеров, поэтому нам как раз помогает п.1)
  3. Далее мы считаем, что клетки, которые алгоритм отнес к одному кластеру являются одним клеточным типом (если кластеров больше, чем типов клеток, то возможно, что 2 разных кластера являются одним типом)
  4. После чего по интенсивности экспрессии поверхностных маркеров мы присваиваем кластеру клеточный тип

Давайте для начала напишем удобную функцию для визуализации наших многомерных данных в пространстве размерностью 2, делать мы это будем при помощи алгоритма t-SNE.

Поиграться с красивой визуализацией можно тут.

In [ ]:
def plot_tsne(data, n_iter=1000, 
              perplexity=40, color=None):
    """
    Функция для отрисовки результатов работы t-SNE
        :param data: таблица объекты х признаки
        :param n_iter: число итераций градиентного спуска,
            может быть полезно увеличить, чтобы получить результаты получше
        :param perplexity: 
        :param color: np.array[n_samples] с переменной,
            которой мы хотим покрасить наши наблюдения
        :return tsne_emb: np.array[n_samples, 2] - результаты работы t-SNE
    """
    
    # Сначала сделаем PCA, так как это хорошее начальное приближение для t-SNE
    # позволит алгоритму быстрее сойтись
    pca = PCA().fit(data)
    pca_embedding = pca.transform(data)
    
    # Запустим t-SNE, он выдаст нам для каждого объекта по 2 числа, 
    # которые мы сможем нарисовать
    tnse = TSNE(n_components=2, init=pca_embedding[:, :2], n_jobs=-1,
                n_iter=n_iter, perplexity=perplexity)
    tsne_embedding = tnse.fit_transform(pca_embedding)
    
    sns.scatterplot(x=tsne_embedding[:, 0],
                    y=tsne_embedding[:, 1],
                    hue=color, palette=sns.color_palette("RdPu")[1:])
    plt.xlabel("$TSNE_1$", fontsize=18)
    plt.ylabel("$TSNE_2$", fontsize=18)
    
    # Вернем также результаты t-SNE, так как, если потом захотим перестроить картинку,
    # в таком случае нам не придется ждать заново, просто нарисуем новую с готовыми данными
    
    return tsne_embedding
In [ ]:
tsne_res = plot_tsne(fc_data_cleared)

Кластеризуйте ваши данные:

  1. Попробуйте методы кластеризации из тех, что мы прошли
  2. Выберите лучший на основании метрики silhouette_score (попробуйте также подобрать гиперпараметры)
  3. Присвойте каждому наблюдению метку класса и нарисуйте график t-SNE, покрасив точки метками кластера
In [ ]:
distances = ["euclidean", "manhattan", "canberra", "cosine"]
linkages = ["average", "single", "complete"]
metrics = ["silhouette"]

algos_selection = algorithm_selection(fc_data_cleared, distances, linkages, metrics, 4)
algos_selection
Out[ ]:
distance linkage silhouette
0 euclidean average 0.585574
1 euclidean single 0.142213
2 euclidean complete 0.614336
3 manhattan average 0.610988
4 manhattan single 0.119002
5 manhattan complete 0.594853
6 canberra average 0.178099
7 canberra single 0.186270
8 canberra complete 0.110311
9 cosine average 0.848678
10 cosine single 0.037162
11 cosine complete 0.838835
In [ ]:
algos_selection.loc[[algos_selection[["silhouette"]].idxmax()[0]]]
Out[ ]:
distance linkage silhouette
9 cosine average 0.848678

Let's try the best parameters combination. Actually the result looks good.

In [ ]:
cluster_attempt_1 = AgglomerativeClustering(n_clusters=4, affinity='cosine', linkage='average').fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_1)
;
Out[ ]:
''

What about the worst parameters combination? Looks like a failed clustering. Viva la hyperparameters search!

In [ ]:
algos_selection.loc[[algos_selection[["silhouette"]].idxmin()[0]]]
Out[ ]:
distance linkage silhouette
10 cosine single 0.037162
In [ ]:
cluster_attempt_2 = AgglomerativeClustering(n_clusters=4, affinity='cosine', linkage='single').fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_2)
;
Out[ ]:
''

Удалось ли вам получить ситуацию, где отдельные группы точек покрашены в один цвет?

Actually, thanks to our function for hyperparameters search if the parameters with best silhouette score are chosen, everything is ok.

Okay then I found more algorithms for clustering and wanted to try several; for some of them searched for the best hyperparameters.

KMeans with hyperparameters selection. It seems like algorithm does not impact the silhouette score.

In [ ]:
kmeans_algs = ["lloyd", "elkan", "auto", "full"]
kmeans_alg_selection = pd.DataFrame(columns=["algorithm", "silhouette"])
for alg in kmeans_algs:
    pred = KMeans(n_clusters=4, algorithm=alg).fit_predict(fc_data_cleared)
    kmeans_sil = silhouette_score(fc_data_cleared, pred)
    new_row = list((alg, kmeans_sil))
    kmeans_alg_selection = kmeans_alg_selection.append(pd.Series(new_row, index=kmeans_alg_selection.columns), ignore_index=True)
kmeans_alg_selection
Out[ ]:
algorithm silhouette
0 lloyd 0.613175
1 elkan 0.613175
2 auto 0.613175
3 full 0.613175
In [ ]:
cluster_attempt_3 = KMeans(n_clusters=4).fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_3)
kmeans_sil = silhouette_score(fc_data_cleared, cluster_attempt_3)
print(f'silhouette score: {kmeans_sil}')
;
silhouette score: 0.6131747085038285
Out[ ]:
''

Birch model. Seemingly good, but silhouette score is not much higher than it was for KMeans.

In [ ]:
cluster_attempt_4 = Birch(n_clusters=4).fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_4)
birch_sil = silhouette_score(fc_data_cleared, cluster_attempt_4)
print(f'silhouette score: {birch_sil}')
;
silhouette score: 0.6143232559027585
Out[ ]:
''

DBSCAN with hyperparameter selection.

In [ ]:
epses = [0.5, 0.7, 1, 1.5, 2, 2.5, 2.6, 2.7, 2.8, 3]
minsamples = [5, 10, 15, 20, 30]
dbscan_alg_selection = pd.DataFrame(columns=["eps", "min_samples", "silhouette"])
for e in epses:
    for m in minsamples:
        pred = DBSCAN(eps=e, min_samples=m).fit_predict(fc_data_cleared)
        dbscan_sil = silhouette_score(fc_data_cleared, pred)
        new_row = list((e, m, dbscan_sil))
        dbscan_alg_selection = dbscan_alg_selection.append(pd.Series(new_row, index=dbscan_alg_selection.columns), ignore_index=True)
dbscan_alg_selection
Out[ ]:
eps min_samples silhouette
0 0.5 5.0 -0.472569
1 0.5 10.0 -0.391862
2 0.5 15.0 -0.324688
3 0.5 20.0 -0.373553
4 0.5 30.0 -0.377765
5 0.7 5.0 -0.014553
6 0.7 10.0 0.115968
7 0.7 15.0 0.136115
8 0.7 20.0 0.170604
9 0.7 30.0 0.146936
10 1.0 5.0 0.202807
11 1.0 10.0 0.477838
12 1.0 15.0 0.487598
13 1.0 20.0 0.470606
14 1.0 30.0 0.280750
15 1.5 5.0 0.525658
16 1.5 10.0 0.595351
17 1.5 15.0 0.592083
18 1.5 20.0 0.589563
19 1.5 30.0 0.583871
20 2.0 5.0 0.574684
21 2.0 10.0 0.565001
22 2.0 15.0 0.599053
23 2.0 20.0 0.598886
24 2.0 30.0 0.598446
25 2.5 5.0 0.608165
26 2.5 10.0 0.609505
27 2.5 15.0 0.608230
28 2.5 20.0 0.604067
29 2.5 30.0 0.605692
30 2.6 5.0 0.504360
31 2.6 10.0 0.521073
32 2.6 15.0 0.521698
33 2.6 20.0 0.503953
34 2.6 30.0 0.504718
35 2.7 5.0 0.470029
36 2.7 10.0 0.500715
37 2.7 15.0 0.517098
38 2.7 20.0 0.515398
39 2.7 30.0 0.498153
40 2.8 5.0 0.481947
41 2.8 10.0 0.510145
42 2.8 15.0 0.516912
43 2.8 20.0 0.520732
44 2.8 30.0 0.505539
45 3.0 5.0 0.266199
46 3.0 10.0 0.266199
47 3.0 15.0 0.266199
48 3.0 20.0 0.442567
49 3.0 30.0 0.478797
In [ ]:
dbscan_alg_selection.loc[[dbscan_alg_selection[["silhouette"]].idxmax()[0]]]
Out[ ]:
eps min_samples silhouette
26 2.5 10.0 0.609505
In [ ]:
cluster_attempt_5 = DBSCAN(eps=2.5, min_samples=10).fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_5)
dbscan_sil = silhouette_score(fc_data_cleared, cluster_attempt_5)
print(f'silhouette score: {dbscan_sil}')
;
silhouette score: 0.6095054495899878
Out[ ]:
''

BisectingKMeans.

In [ ]:
cluster_attempt_6 = BisectingKMeans(n_clusters=4).fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_6)
bisect_sil = silhouette_score(fc_data_cleared, cluster_attempt_6)
print(f'silhouette score: {bisect_sil}')
;
silhouette score: 0.610081988930565
Out[ ]:
''

MiniBatchKMeans.

In [ ]:
cluster_attempt_7 = MiniBatchKMeans(n_clusters=4).fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_7)
minibatch_sil = silhouette_score(fc_data_cleared, cluster_attempt_7)
print(f'silhouette score: {minibatch_sil}')
;
silhouette score: 0.613694891921406
Out[ ]:
''

SpectralClustering with hyperparameters selection.

(runs very slowly due to lots of hyperparameter combinations)

In [ ]:
es = ["arpack", "lobpcg", "amg"]
aff = ['cosine', 'chi2', 'linear', 'rbf', 'nearest_neighbors']
nneighb = [5, 10, 15, 20, 30]
als = ["kmeans", "discretize", "cluster_qr"]
specclus_alg_selection = pd.DataFrame(columns=["eigen_solver", "affinity", "n_neighbors", "assign_labels", "silhouette"])
for e in es:
    for a in aff:
        for n in nneighb:
            for al in als:
                try:
                    pred = SpectralClustering(n_clusters=4, eigen_solver=e, 
                                            affinity=a, n_neighbors=n, assign_labels=al).fit_predict(fc_data_cleared)
                    specclus_sil = silhouette_score(fc_data_cleared, pred)
                    new_row = list((e, a, n, al, specclus_sil))
                    specclus_alg_selection = specclus_alg_selection.append(pd.Series(new_row, index=specclus_alg_selection.columns), ignore_index=True)
                except ValueError:
                    continue
specclus_alg_selection
Out[ ]:
eigen_solver affinity n_neighbors assign_labels silhouette
0 arpack cosine 5 kmeans 0.614579
1 arpack cosine 5 discretize 0.579351
2 arpack cosine 5 cluster_qr 0.579911
3 arpack cosine 10 kmeans 0.614579
4 arpack cosine 10 discretize 0.579351
... ... ... ... ... ...
115 lobpcg nearest_neighbors 20 discretize 0.614323
116 lobpcg nearest_neighbors 20 cluster_qr 0.614323
117 lobpcg nearest_neighbors 30 kmeans 0.614323
118 lobpcg nearest_neighbors 30 discretize 0.614423
119 lobpcg nearest_neighbors 30 cluster_qr 0.614423

120 rows × 5 columns

In [ ]:
specclus_alg_selection.loc[[specclus_alg_selection[["silhouette"]].idxmax()[0]]]
Out[ ]:
eigen_solver affinity n_neighbors assign_labels silhouette
15 arpack linear 5 kmeans 0.614694
In [ ]:
cluster_attempt_8 = SpectralClustering(n_clusters=4, eigen_solver="arpack", affinity="linear", 
                                       n_neighbors=5, assign_labels="kmeans").fit_predict(fc_data_cleared)
plot_tsne(fc_data_cleared, color=cluster_attempt_8)
specclus_sil = silhouette_score(fc_data_cleared, cluster_attempt_8)
print(f'silhouette score: {specclus_sil}')
;
silhouette score: 0.6146941826291717
Out[ ]:
''

3.3. Аннотация клеточных типов¶

4 балла

Теперь когда мы подобрали хороший алгоритм кластеризации, можно аннотировать наши клетки. Для этого мы нарисуем t-SNE и покрасим точки в зависимости от интенсивности экспрессии поверхностных маркеров. В датасете присутствуют следующие типы клеток:

  • B_cells
  • T_cells
  • Monocytes
  • Other cells

Вам нужно будет выяснить, какие клетки экспрессируют определенные маркеры и присвоить каждому кластеру один из типов клеток.

Для начала нарисуем все это безобразие:

In [ ]:
# Результаты t-SNE уже есть в переменной tsne_res
fig, axes = plt.subplots(3, 3, figsize=(20, 20))
for col, ax in zip(fc_data_cleared.columns, axes.ravel()):
    scatter = ax.scatter(tsne_res[:, 0], tsne_res[:, 1], 
                         c=fc_data_cleared[col], cmap="RdPu")
    fig.colorbar(scatter, ax=ax)
    ax.set_title(col)
    ax.grid(False)
    ax.axis("off")

Дальше дело за вами, нужно определить клеточный тип для каждого кластера и записать их как значения в словаре:

lineage markers for T cells (CD3), B cells (CD19) and monocytes (CD14)

In [ ]:
cell_type_cluster_map = {0: "T_cell", 
                         1: "B_cell", 
                         2: "monocyte", 
                         3: "other"}
In [ ]:
clust = AgglomerativeClustering(n_clusters=4, affinity='cosine', linkage='average').fit(fc_data_cleared)
In [ ]:
labeled_fc_data = fc_data_cleared.assign(Population=clust.labels_)
labeled_fc_data["Population"] = labeled_fc_data["Population"].map(cell_type_cluster_map)

Посчитайте, как распределены типы клеток:

In [ ]:
labeled_fc_data.head()
Out[ ]:
FSC-A SSC-A HLA-DR CD3 CD45 CD14 CD56 CD19 CD11c Population
0 0.722199 0.154232 0.494449 0.512853 4.057328 1.019255 4.007273 -1.406912 0.384406 other
1 0.751073 0.114618 4.671924 0.822042 4.933149 0.159388 -0.216407 6.034845 0.577039 monocyte
2 1.108188 0.150011 -0.422857 4.006353 5.616384 0.519223 0.783077 -0.318668 -0.493832 T_cell
3 0.553310 0.119929 0.345968 3.332310 5.260935 1.481229 0.264846 0.694883 3.416960 T_cell
4 0.774884 0.076984 4.984963 0.885501 4.562367 0.217117 0.086279 5.563345 0.471591 monocyte
In [ ]:
labeled_fc_data["Population"].value_counts() # class disbalance detected
Out[ ]:
T_cell      6222
monocyte    2035
B_cell      1100
other        632
Name: Population, dtype: int64

Сохраните ваши результаты в csv файл, мы сравним их с правильными ответами по метрике f1_score(average="macro"), которая хорошо подходит, даже если классы не сбалансированы.

In [ ]:
labeled_fc_data.to_csv("../data/labeled_cells.csv")

3.4. Атоматическая аннотация¶

дополнительные 5 баллов

Это задание является дополнительным. Вам нужно построить классификатор клеточных типов на основе размеченных вами данных. Пока что мы прошли не так много разных классификаторов, но уже есть из чего выбирать.

Оцените качество полученной модели при помощи f1_score(average="macro"). Попробуйте его улучшить, каким-то образом поработав с дисбалансом классов (помощь).

First attempt. Without class disbalance correction.

In [ ]:
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.utils import resample
from sklearn.preprocessing import StandardScaler
In [ ]:
# train-test split

X = labeled_fc_data.drop("Population", axis=1)
y = labeled_fc_data.Population
test_proportion = 0.25
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_proportion, random_state=SEED)
In [ ]:
# check that proportion of each class in train and test remains at the setup level (we remember that train_test_split is a random thing)
# see class disbalance, have to fix it for train

frame = {'train': y_train.value_counts(), 'test': y_test.value_counts()}
result = pd.DataFrame(frame)
result["proportion"] = result.test / (result.train + result.test)
result
Out[ ]:
train test proportion
T_cell 4696 1526 0.245259
monocyte 1513 522 0.256511
B_cell 812 288 0.261818
other 470 162 0.256329

The major class in train is T-cell (4696 cells). We have to upsample other three classes up to this number.

In [ ]:
train_overall = pd.concat([X_train, y_train], axis=1)

train_all_mon = train_overall[train_overall["Population"] == "monocyte"]
train_all_B = train_overall[train_overall["Population"] == "B_cell"]
train_all_other = train_overall[train_overall["Population"] == "other"]
train_all_T = train_overall[train_overall["Population"] == "T_cell"]

reshape_val = train_all_T.shape[0]

train_all_mon_up = resample(train_all_mon, random_state=SEED, n_samples=reshape_val,replace=True)
train_all_B_up = resample(train_all_B, random_state=SEED, n_samples=reshape_val,replace=True)
train_all_other_up = resample(train_all_other, random_state=SEED, n_samples=reshape_val,replace=True)

train_balanced = pd.concat([train_all_B_up, train_all_T, train_all_mon_up, train_all_other_up])
train_balanced.Population.value_counts()
Out[ ]:
B_cell      4696
T_cell      4696
monocyte    4696
other       4696
Name: Population, dtype: int64
In [ ]:
X_train, y_train = train_balanced.drop("Population", axis=1), train_balanced.Population

Now I want to standardize data.

In [ ]:
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
pd.DataFrame(X_train_scaled).describe()
Out[ ]:
0 1 2 3 4 5 6 7 8
count 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04 1.878400e+04
mean -1.331511e-16 -3.631394e-17 1.210465e-16 -6.052323e-17 -6.052323e-17 -1.210465e-17 -4.841859e-17 -1.210465e-17 1.210465e-16
std 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00 1.000027e+00
min -2.833427e+00 -1.198775e+00 -3.230950e+00 -1.958418e+00 -7.761042e+00 -1.203148e+00 -1.700852e+00 -1.639406e+00 -2.074041e+00
25% -6.274094e-01 -6.331310e-01 -9.057285e-01 -7.428025e-01 -5.894646e-01 -6.454102e-01 -7.230285e-01 -6.542683e-01 -7.749901e-01
50% -2.485926e-01 -4.030355e-01 -7.330673e-02 -2.865904e-01 1.172072e-02 -4.689578e-01 -4.941840e-01 -4.531557e-01 -2.333854e-01
75% 3.941363e-01 1.139127e-01 9.313643e-01 1.009500e+00 6.915309e-01 9.250426e-02 9.817381e-01 5.582255e-01 6.737365e-01
max 9.666798e+00 5.631934e+00 2.068627e+00 3.481183e+00 3.630924e+00 2.505861e+00 3.158167e+00 2.321100e+00 2.553138e+00
In [ ]:
def find_best_k_accuracy(X_train, y_train, X_test, y_test, params, beau=False):
    """
    Choose the best k for KKNClassifier
    Arguments:
    X_train, np array (num_train_samples, num_features) - train data
    y_train, np array (num_train_samples) - train labels
    X_test, np array (num_test_samples, num_features) - test data
    y_test, np array (num_test_samples) - test labels
    params, list of hyperparameters for KNN, here it is list of k values
    metric, function for metric calculation
    Returns:
    train_metrics the list of metric values on train data set for each k in params
    test_metrics the list of metric values on test data set for each k in params
    """
    
    """
    YOUR CODE IS HERE
    """

    result_train = list()
    result_test = list()

    for neighbour_number in params:
        knn_classifier_search = KNeighborsClassifier(n_neighbors=neighbour_number, weights="distance")
        knn_classifier_search.fit(X_train, y_train)

        prediction_search_train = knn_classifier_search.predict(X_train)
        prediction_search_test = knn_classifier_search.predict(X_test)

        metric_train = accuracy_score(y_train, prediction_search_train)
        metric_test = accuracy_score(y_test, prediction_search_test)

        if beau:
            print(f'K number: {neighbour_number}')
            print(f'Train dataset metric: {metric_train}')
            print(f'Test dataset metric: {metric_test}')

        result_train.append(metric_train)
        result_test.append(metric_test)

    return result_train, result_test
In [ ]:
params = list(range(1, 31))
train_metrics, test_metrics = find_best_k_accuracy(X_train_scaled, y_train, X_test, y_test, params, beau=False)

Optimal neighbour number: 2, 11, 12.

In [ ]:
plt.title("Accuracy for optimal neighbour number search in multiclass KNN classification model", size=15, weight="bold")
plt.plot(params, train_metrics, label="train", color=sns.color_palette("RdPu")[5])
plt.plot(params, test_metrics, label="test", color=sns.color_palette("RdPu")[2])
plt.xlabel("K in KNN", size=13, weight="bold")
plt.ylabel("Accuracy", size=13, weight="bold")
plt.axvline(x=2, color=sns.color_palette("RdPu")[4], linestyle='-.', linewidth=1.2)
plt.text(x=1, y=0.99, s="N=2", fontdict={'fontsize': 10, 'fontweight': 'bold', 'color': sns.color_palette("RdPu")[4]})
plt.axvline(x=11, color=sns.color_palette("RdPu")[4], linestyle='-.', linewidth=1.2)
plt.text(x=9.8, y=0.99, s="N=11", fontdict={'fontsize': 10, 'fontweight': 'bold', 'color': sns.color_palette("RdPu")[4]})
plt.axhline(y=max(test_metrics), color=sns.color_palette("RdPu")[4], linestyle='-.', linewidth=1.2)
plt.text(x=-1.8, y=max(test_metrics), s=max(test_metrics).round(3), fontdict={'fontsize': 10, 'fontweight': 'bold', 'color': sns.color_palette("RdPu")[4]})
plt.legend(loc='upper left')
;
Out[ ]:
''
In [ ]:
mdl = KNeighborsClassifier(n_neighbors=11, weights="distance")
mdl.fit(X_train_scaled, y_train)
len(y_test) - sum(mdl.predict(X_test) == y_test)
Out[ ]:
166

Задание 4. Любимое место¶

дополнительно 0.5 балла

Расскажите о вашем любимом месте, куда бы вы хотели съездить этой весной, если вы уже там бывали, то прикрепите фотографию)

Actually my plans might be pretty spontaneous this spring. Hope that I'll make it to Saint Petersburg to celebrate the BI's B-day. Or maybe Yerevan of Istanbul in blossom of the first flowers at the end of May? Tons of dreams, tons of plans, tons of work at the same time :)

Anyway, I attach a picture for each place.

In [ ]:
from IPython.display import Image
In [ ]:
Image("spb.jpg") # one of the most beautiful locations in St Petersburg
Out[ ]:
In [ ]:
Image("ist.jpg") # a piece of good Istanbul memories
Out[ ]:
In [ ]:
Image("ere.jpg") # tiles from Yerevan (wish my github commits calendar would look like this)
Out[ ]:

Therapy time¶

Напишите здесь ваши впечатления о задании: было ли интересно, было ли слишком легко или наоборот сложно и тд. Также сюда можно написать свои идеи по улучшению заданий, а также предложить данные, на основе которых вы бы хотели построить следующие дз.

Ваши мысли:

In [ ]: